Diffusion of small light particles in a solvent of large massive molecules 
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We study diffusion of small light particles in a solvent which consists of large heavy particles. 
The intermolecular interactions are chosen to approximately mimic a water-sucrose (or water- 
polysaccharide) mixture. Both computer simulation and mode coupling theoretical (MCT) cal- 
culations have been performed for a solvent-to-solute size ratio five and for a large variation of the 
mass ratio, keeping the mass of the solute fixed. Even in the limit of large mass ratio the solute 
motion is found to remain surprisingly coupled to the solvent dynamics. Interestingly, at intermedi- 
ate values of the mass ratio, the self-intermediate scattering function of the solute, F s (k,t) (where 
k is the wavenumber and t the time), develops a stretching at long time which could be fitted to a 
stretched exponential function with a k-dependent exponent, /3. For very large mass ratio, we find 
the existence of two stretched exponentials separated by a power law type plateau. The analysis of 
the trajectory shows the coexistence of both hopping and continuous motions for both the solute 
and the solvent particles. It is found that for mass ratio five, the MCT calculations of the self- 
diffusion underestimates the simulated value by about 20%, which appears to be reasonable because 
the conventional form of MCT does not include the hopping mode. However, for larger mass ratio, 
MCT appears to breakdown more severely. The breakdown of the MCT for large mass ratio can be 
connected to a similar breakdown near the glass transition. 



I. INTRODUCTION 

The issue of diffusion of small light particles in a sol- 
vent composed of larger and heavier particles is uncon- 
ventional because the role of the solvent in the solute 
diffusion is different from the case where the sizes are 
comparable. There are two limits that can be identified 
for such systems. One limit is the well studied Lorentz 
gas system, which consists of a single point particle mov- 
ing in a triangular array of immobile disk scatters. Here 
the motion of the point particle can be modelled by ran- 
dom walk between traps. 1-3 The other limit is where the 
size of the solute particle is still smaller than that of the 
solvent molecules but it has a finite size (that is, not 
a point) and while the solvent is slow (compared to the 
solute particles) but not completely immobile. In the lat- 
ter case, the translational diffusion of the solute is often 
attempted to describe by the well known hydrodynamic 
Stokes-Einstein (SE) relation given by, 4,5 



D = 



k B T 
CtjR 



(1) 



where ks is the Boltzmann constant, T is the absolute 
temperature, C is a numerical constant determined by 
the hydrodynamic boundary condition, rj is the shear vis- 
cosity of the solvent and R is the radius of the diffusing 



particle. Validity of Eq. 1 for small solutes is, of course, 
questionable. 4-6 

There have been many experimental, 6,7 computer 
simulation 8-10 and theoretical 11,12 studies of diffusion 
of small solute particles in a solvent composed of larger 
particles. All these studies show that the SE relation 
significantly underestimates the diffusion coefficient. To 
explain the enhanced diffusion, sometimes an empirical 
modification of the SE relation is used. 6 ' 7 It is considered 
that D oc ?y -Q , where a ~ 2/3. This fractional viscos- 
ity dependence is referred to as the microviscosity effect 
which implies that the viscosity around the small solute 
is rather different from that of the bulk viscosity. The en- 
hanced diffusion value has also been explained in terms 
of effective hydrodynamic radius. ' 13 

The earlier mode coupling theoretical (MCT) 
studies 11 ' 12 of diffusion of smaller solutes in a solvent 
composed of larger size molecules attributed the en- 
hanced diffusion to the decoupling of the solute motion 
from the structural relaxation of the solvent. The MCT 
studies suggest that this decoupling of the solute motion 
from the structural relaxation of the solvent can lead 
to the fractional viscosity dependence often observed in 
supercooled liquids. 

However, there have been no systematic studies of the 
effects of the variation of size and mass of the solute- 
solvent system. In this article we have explored the dif- 
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fusional mechanism of the isolated small particles (so- 
lute) in a liquid composed of larger particles (solvent), 
both analytically and numerically. The study is per- 
formed by keeping the solvent-to-solute size ratio (Sr) 
fixed at five, but varying the mass of the solvent over 
a large range, by keeping the mass of the solute fixed. 
That is, the mass ratio Mr (solvent mass/solute mass) 
is progressively raised to higher values. This system is 
expected to mimic some aspects of the water-sucrose or 
water-polysaccharide mixtures. 14 

The trajectories of the solute and the solvent show 
coexistence of both hopping and continuous motions. 
As the solvent mass is increased, the self-intermediate 
scattering function of the solute develops an interesting 
stretching in the long time. For larger mass ratio we see 
the existence of two stretched exponentials separated by 
a power law type plateau. 

The mode coupling theory calculation of the self- 
diffusion coefficient of solute particles performed in the 
limit of small mass ratio (of 5) is found to be in qualita- 
tive agreement with the simulated diffusion - MCT un- 
derestimates the diffusion by about 20%. Thus, although 
the MCT underestimates the diffusion, the agreement is 
satisfactory in light of the contribution from the hop- 
ping mode to diffusion which MCT does not explicitly 
take into account. However, the deviation from the sim- 
ulated value increases with increase in mass ratio (which 
is equivalent to the increase of the mass of the solvent 
particles). In the limit of large mass ratio, MCT breaks 
down. The binary contribution to the total friction is 
found to decrease as one increases the mass of the sol- 
vent. In addition, due to the development of stretching 
in the self-intermediate scattering function of the solute 
and the inherent slow solvent dynamics, there remains a 
strong coupling of the solute motion to the solvent den- 
sity mode. This enhanced coupling at larger mass ratio 
came as a surprise to us. 

In the limit of very large mass ratio, the motion of the 
light solute particle resembles to that of its motion in an 
almost frozen disorder system like near the glass transi- 
tion temperature. Thus, the breakdown of MCT in the 
limit of large mass ratio could be connected to its fail- 
ure near the glass transition temperature. Of course, one 
should note that in the limit of mass of the solvent goes 
to infinity, the basic assumption of conventional MCT 
breaks down. 

It is widely believed that in a deeply supercooled liquid 
close to its glass transition temperature (T g ), the hopping 
mode is the dominant mode in the system which controls 
the mass transport and the stress relaxation. Recently, 
a computer simulation study of a deeply supercooled bi- 
nary mixture 15 has shown the evidence of an intimate 
connection between the anisotropy in local stress and the 
particle hopping. It was shown that the local anisotropy 
in the stress is responsible for the particle hopping in a 
particular direction. Furthermore, it was suggested that 
the local frustration present in the system (which is more 
in a binary mixture with components of different sizes) 



could cause the local anisotropy in the stress which in 
turn acts as a driving force for hopping. However, in the 
present study, the density (or the pressure) of the system 
is not as high as that of a deeply supercooled system. The 
relaxation of the stress is found to occur much faster and 
it relaxes almost completely within our simulation time 
window even for the largest mass ratio. Consequently, 
the microscopic origin of particle hopping here could be 
different than that for a deeply supercooled liquid. 

The layout of the rest of the paper is as follows. Sec- 
tion II deals with the system and simulation details. The 
simulation results and discussions are given in the next 
Sec. Ill and the mode coupling theoretical analysis is 
presented in Sec. IV. In Sec. V, we discuss the possi- 
ble effect of heterogeneity probed by the solute on the 
self-dynamic intermediate scattering function of the so- 
lute. Finally, in Sec. VI we present the conclusions of 
the study. 



II. SYSTEM AND SIMULATION DETAILS 

We performed a series of equilibrium isothermal- 
isobaric (N P T) ensemble molecular dynamics (MD) 
simulation of binary mixtures in three dimensions for an 
infinitesimal small value of the mole fraction of one of 
the species. The binary system studied here contains a 
total of TV = 500 particles consisting of two species of 
particles, with N\ = 490 and 7V2 = 10 number of parti- 
cles. Hereafter, we refer the indices I and 2, respectively, 
for the solvent and solute particles. Thus, the mixture 
under study contains 2% of solute particles. The inter- 
action between any two particles is modeled by means of 
shifted force Lennard- Jones (LJ) pair potential, 16 where 
the standard LJ is given by 
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where i and j denote two different particles (1 and 2). 
In our model system, the potential parameters are as 
follows: en = 1.0, an = 1.0, e 2 2 = 1.0, 022 = 0.2, 
612 = 2.0 (enhanced attraction), and u\i = 0.6. The 
mass of the solute particles is chosen to be m 2 = 0.2 
where the solvent (species 1) mass mi is increasingly var- 
ied and four different values are chosen 1, 5, 10 and 50. 
Thus, in this study we examined four different solvent- 
to-solute mass ratio, Mr — mi/mz — 5, 25, 50 and 250 
for a fixed solvent-to-solute size ratio, Sr = owloii = 5. 
Note that in the model system being studied the solute- 
solvent interaction (ei 2 ) is much stronger than both of 
their respective pure counterparts. In order to lower the 
computational burden the potential has been truncated 
with a cutoff radius of 2.5<7n. All the quantities in this 
study are given in reduced units, that is, length in units 
of (Tn, temperature T in units of eu/ks, pressure P in 
units of en/<7ii, and the mass in the unit of m, which can 
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be assumed as argon (Ar) mass unit. The corresponding 
microscopic time scale is r = y/maf 1 /eu. 

All simulations in the NPT ensemble were performed 
using the Nose- Hoover- Andersen method, 17 where the 
external reduced temperature is held fixed at T* = 0.8. 
The external reduced pressure has been kept fixed at 
P* = 6.0. The reduced average density p* of the sys- 
tem corresponding to this thermodynamic state point is 
0.989 for all the mass ratios being studied. Through- 
out the course of the simulations, the barostat and the 
system's degrees of freedom are coupled to an indepen- 
dent Nose-Hoover chain 18 (NHC) of thermostats, each of 
length 5. The extended system equations of motion are 
integrated using the reversible integrator method 19 with 
a small time step of 0.0002. The higher order multiple 
time step method 20 has been employed in the NHC evo- 
lution operator which lead to stable energy conservation 
for non-Hamiltonian dynamical systems. 21 The extended 
system time scale parameter used in the calculations was 
taken to be 0.9274 for both the barostat and thermostats. 

The systems were equilibrated for 10 6 time steps and 
simulations were carried out for another 2.0 x 10 6 pro- 
duction steps following equilibration, during which the 
quantities of interest are calculated. The dynamic quan- 
tities are averaged over three such independent runs for 
better improvement of the statistics. We have also calcu- 
lated the partial radial distribution functions (ffn(r) and 
912(f)) to make sure that there is no crystallization. 



III. SIMULATION RESULTS AND DISCUSSION 

In figure 1 we show typical solute trajectories for four 
different solvent-to-solute mass ratio, Mr (= mijmi-, vni 
is the mass of the solvent particles). The trajectories 
reveal interesting dependence on Mr. At the value of 
Mr equal to five, the solute trajectory is mostly con- 
tinuous with occasional hops. As the mass ratio Mr is 
increased, the solute motion gets more trapped and its 
motion tends to become discontinuous where displace- 
ments occur mostly by hopping. This is because with 
increase in the solvent mass the time scale of motion of 
the solvent particles become increasingly slower. Thus, 
the solute gets caged by the solvent particles and keeps 
rattling between a cage till one solvent particle moves 
considerably to disperse the solute trajectory (see trajec- 
tory for Mr = 250). Thus, there is a remarkable change 
in solute's motion in going from Mr = 5 (figure la) to 
Mr = 250 (figure Id). 

In figure 2 we plot the solvent trajectories for differ- 
ent solvent-to-solute mass ratio, Mr. We find that for 
all values of Mr, there is a coexistence of hopping and 
continuous motion of the solvent molecules. At higher 
solvent mass, as expected, the magnitude of displace- 
ment becomes less and hopping becomes less frequent, 
but, surprisingly, the jump motion still persists. 

Figure 3 displays the decay behavior of the self- 



intermediate scattering function (F s (k,t)) of the solute 
for different mass ratio Mr, at reduced wavenumbcr 
k* = ka\\ ~ 27r. The plot shows that F s (k,t) begins 
to stretch more for higher solvent mass. This stretching 
of F s (k,t) is kind of novel and we have examined it in 
detail. 

After the initial Gaussian decay, F s (k,t), for smaller 
values of Mr, can be fitted to a single stretched expo- 
nential where the exponent (5 ~ 0.6. However, for higher 
mass of the solvent, F s (k,t) can be fitted only to a sum 
of two stretched exponentials. 

The behavior of F s (k,t) of the solvent for all the sol- 
vent masses studied is shown in figure 4. The plot shows 
that (as expected) the time scale of relaxation of F s (k, t) 
becomes longer as the solvent mass is increased. How- 
ever, the self-intermediate scattering function of the sol- 
vent does not display any stretching at long times, even 
for the largest mass ratio considered. The decay can be 
fitted by sum of a short-time Gaussian and a long-time 
exponential function. 

The reason that F s (k,t) of the solute shows such 
stretching but that of the solvent does not, can be ex- 
plained as follows. Due to the small size and the lighter 
weight of the solute, the time scale of motion of the solute 
particles is much shorter compared to that of the solvent 
particles. Consequently, the solute motion probes more 
heterogeneity during the time scale of decay of F s (k, t) of 
the solute. This heterogeneity probed by the solute in- 
creases as the solvent mass is increased. Since the solvent 
motion is much slower, it probes enough configurations 
during the time scale of decay of its F s (k, t). 

In order to quantify the degree of heterogeneity probed 
by the solute, we have plotted the non-Gaussian param- 
eter a2(t), 22 for the solute, in figure 5. Clearly, the het- 
erogeneity probed by the solute (quantified by the peak 
height of a 2 (t)) increases as the solvent mass is increased. 
On the other hand, a 2 (t) of the solvent shows no such 
increase in the peak height of a 2 (t) which remains small 
and unaltered, although the position of the peak shifts to 
longer time as the mass of the solvent is increased. We 
have discussed this analysis in more detail in section V. 

The role of the local heterogeneity can be further ex- 
plored by calculating F s (k, t) at wavenumber correspond- 
ing to solute-solvent average separation. This corre- 
sponds to k* = lit 1 011, where CT12 = \(o~\-\-o~ 2 ). In figure 
6 we plot the self dynamic structure factor of the so- 
lute for different mass ratio Mr, at reduced wavenumbcr 
k* ~ 27r/<7i2. This is primarily the wavenumber probed 
by the solute. At this wavenumbcr one observes more 
stretching of F s (k, t) at longer times, for all the mass ra- 
tios. This is in agreement with the above argument that 
since the time window probed by the solute is smaller at 
higher k, it probes even larger heterogeneity. 

The emergence of the plateau between the two 
stretched exponentials in F s (k,t) of the solute (figure 6), 
can be attributed to the separation of time scale of the 
binary collision and the solvent density mode contribu- 
tion to the particle motion. This separation of time scale 
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increases as the solvent mass is increased. As the decay 
after the plateau is mainly due to the density mode con- 
tribution, the plateau becomes more prominent as the 
mass of the solvent is increased. Table I clearly shows 
the separation of time scales where the value of the time 
constants and the exponents obtained from the two dif- 
ferent stretched exponential fits to F s (k, t) are presented 
for different mass ratio Mr. 

The F s (k,t) of each of the individual solute particle 
obtained from a single MD run is shown in figure 7, at 
k* ~ 27t/cti2 and for Mr = 250. We find that not only 
all of them have different time scale of relaxation but 
each of them shows considerable stretching (the expo- 
nent (3 ~ 0.6) at longer time. This confirms further that 
each of the solute particle probes heterogeneous structure 
and dynamics of the solvent. 

In table II we present the scaled average (over all the 
solute particles and three independent MD runs) diffu- 
sion value of the solute particles obtained from the slope 
of the mean square displacement (MSD) in the diffusive 
limit, for different mass ratio, Mr. The values of the 
solute diffusion decreases as the mass of the solvent is in- 
creased, as expected. The mass dependence can be fitted 
to a power law as clearly manifested in figure 8 where we 
have plotted In 1/D 2 against In mi/m2, where D2 is the 
self-diffusion of the solute particles. The slope of the line 
is about 0.13. The small value of the exponent is clearly 
an indication of the weak mass dependence of the self- 
diffusion coefficient of small solute particle on the mass 
of the bigger solvent particles. 

Interestingly, it is to be noted that a similar weak 
power law mass dependence was seen in the self-diffusion 
coefficient of a tagged particle on its mass - the expo- 
nent was often found to be around 0.1. Recently, a self- 
consistent mode coupling theory (MCT) analysis success- 
fully explained this weak mass dependence. 23 

In figure 9 we plot the normalized velocity autocorrela- 
tion function C v (t) of the solute particles for the different 
values of Mr. The velocity correlation function shows 
highly interesting features at larger mass ratio. Not 
only does the negative dip becomes larger, there devel- 
ops a second minimum or an extended negative plateau 
which becomes prominent as the mass of the solvent is 
increased. Interestingly, as can be seen from the figure 
the C v (t) shows an oscillatory behavior that persists for 
long period. This is clearly evidence for the 'dynamic 
cage' formation in which the solute particle is seen to 
execute a damped oscillatory motion. Because of the 
increasing effective structural rigidity of the neighbor- 
ing solvent particles as the mass of the solvent increases, 
the motion of the solute particle can be modelled as a 
damped oscillator which is reminiscent of the behavior 
observed in a deeply supercooled liquid near the glass 
transition temperature. 24 

The understanding of the microscopic origin of the de- 
velopment of an increasingly negative dip followed by 
pronounced oscillations at longer times in the velocity au- 
tocorrelation function of a supercooled liquid is a subject 



of much current interest. The novel molecular-dynamics 
simulation study of Kivelson and his coworkers 25 surely 
provide a step forward in this direction. Their simulation 
study had shown that the single particle velocity auto- 
correlation function could be thought of as a sum of the 
local rattling motion relative to the center of mass of the 
neighboring cluster and the motion of the center of mass 
of that cluster. Furthermore, it was observed that the 
rich structure displayed by the velocity autocorrelation 
function at high density arise primarily from the relative 
motion of the tagged particle, that is, the rattling motion 
within the cage formed by the neighboring particles. 



IV. MODE COUPLING THEORY ANALYSIS 

Mode coupling theory remains the only quantitative 
fully microscopic theory for self-diffusion in strongly cor- 
related systems. In this section we present a mode cou- 
pling theory calculation of the solute diffusion for differ- 
ent solvent-to-solute mass ratio Mr. Diffusion coefficient 
of a tagged solute is given by the well known Einstein re- 
lation, 



D 2 = k B T/m 2 ( 2 (z = 0), 



(3) 



where D 2 is the diffusion coefficient of the solute and 
(2(2) is the frequency dependent friction. m 2 is the mass 
of the solute particle. Mode coupling theory provides 
an expression of the frequency dependent friction on an 
isolated solute in a solvent. 

In the normal liquid regime (in the absence of hopping 
transport) it can be given by, 11,12 



1 
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C2W 



ti(z) + R%(z) 



R^(z) 



(4) 



where C 2 ( z ) is the binary part of the friction, i?2i( z ) i s 
the friction due to the coupling of the solute motion to 
the collective density mode of the solvent and R 2 ^{z) is 
the contribution to the diffusion (inverse of friction) from 
the current modes of the solvent. 

For the present system, we have neglected the contri- 
bution from the current term, (z) which is expected 
to be reasonable at high density and low temperature. 
Thus the total frequency dependent friction can be ap- 
proximated as, 



( 2 (z) ~ C 2 B (z) + R p £{z) 



(5) 



The expression for the time dependent binary friction 
( 2 {t), for solute-solvent pair, is given by, 11 ' 12 



C2 (*) = Lo 2 ol2 exp(-t'/T^) 



(6) 



where lo i 2 is now the Einstein frequency of the solute in 
presence of the solvent and is given by, 



w ol2 = 



P 

3m 2 



J dvg 12 {r)V 2 v 12 (r). (7) 
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Here gi 2 (r) is the partial solute-solvent radial distribu- 
tion function. 

In Eq. 6, the relaxation time is determined from the 
second derivative of ( 2 (t) at t = and is given by, 12 ' 23 



^012 



/r 2 = (p/6m 2f i) / dv{V a V f3 v 12 {v))g l2 {r) 

x (V tt V^ 12 (r)) + (l/6p) J [dk/(2n) 3 ] 
x 7d f 2 (k)(5(fc)-l) 7d f 2 (k) (8) 



where summation over repeated indices is implied. \i is 
the reduced mass of the solute-solvent pair. Here S(q) 
is the static structure factor which is obtained from the 
HMSA scheme. 26 The expression for 7 ^(k) is written 
as a combination of the distinct parts of the second mo- 
ments of the longitudinal and transverse current correla- 
tion functions 7^ 12 (k) and 7^12 ( k )) respectively. 

Idi2( k ) = -{p/ m 2) j dr exp(-ik.r),gi 2 (r)V a V /3 wi2(r) 



k a k^ dl2 (k) + (S a ^k a k' 3 ) 7dl2 (k) 



(9) 



wh ere 7di 2 ( k ) = 7^2 ( k ) and 7di 2 ( k ) = 7?f 2 ( k )- 

The expression for R 2 1(t), for solute-solvent pair, can 
be written as, 12 ' 23 

RPP{t) = pksT r [ M/(2 w )S]$.k>) 2 k>*[c 12 {k>)] 2 
m 2 J 

x [F s (k',t)F{k',t) - F s a {k' ,t)F {k' ,t)] (10) 

In Eq. 10, c\ 2 {k) is the two particle (solute-solvent) 
direct correlation in the wavenumber (fc) space which is 
obtained here from the HMSA scheme. 26 The partial ra- 
dial distribution function (gi 2 (r)) required to calculate 
the Einstein frequency (lv i 2 ) and the binary time con- 
stant (t£ ) is obtained from the present simulation study. 
F(k, t) is the intermediate scattering function of the sol- 
vent, and F (k,t) is the inertial part of the intermediate 
scattering function. F s (k, t) is the self-intcrmcdiate scat- 
tering function of the solute and F£(k,i) is the inertial 
part of F s (k, t). 

It should be noted here that the short time dynam- 
ics of the density term used in Eq.10, is different from 
the conventional mode coupling formalism. 27 This pre- 
scription has recently been proposed to explain the weak 
power law mass dependence of the self-diffusion coeffi- 
cient of a tagged particle. 23 The detailed discussion on 
this prescription has been given elsewhere. 12 ' 23 

Since the solvent is much heavier than the solute, 
the decay of solvent dynamical variables are naturally 
much slower than those of the solute. Since the decay of 
F s (k,t) is much faster than F(k,t), in Eq. 10 the con- 
tribution from the product, F s (k,t)F(k,t) mainly gov- 
erned by the time scale of decay of F s (k,t). Thus the 
long time part of the F(k, t) becomes unimportant and 
the viscoelastic expression 12 for F(k, z) (Laplace trans- 
form of F(k,t)) obtained by using the well-known Mori 



continued-fraction expansion, truncating at second order 
would be a reasonably good approximation. The expres- 
sion of F(k, z), can be written as, 11,12 



F(k,z) = 



S(k) 



z + 



(11) 



z + 



+ r k - 



where F(k, t) is obtained by Laplace inversion of F(k, z), 
the dynamic structure factor. Because of the viscoelastic 
approximation, < ui\ > and and also Tk are deter- 
mined by the static pair correlation functions. The static 
pair correlation functions needed are the static structure 
factor S(q) and the partial solvent-solvent radial distribu- 
tion function gn(r). S(q) is obtained by using the HMSA 
scheme 26 and #11 (r) is taken from the present simulation 
study. 

We have used the recently proposed 12 ' 23 generalized 
self-consistent scheme to calculate the friction, (,(z), 
which makes use of the well-known Gaussian approxi- 
mation for F s (k,t), 28 



F s (k,t)= exp( ) 



= exp 



k B T ] 
m 2 



' f drC v {T){t-T) 
Jo 



(12) 



where < Ar 2 (t) > is the mean square displacement 
(MSD) and C v (t) is the time-dependent velocity au- 
tocorrelation function (VACF) of the solute particles. 
The time-dependent VACF is obtained by numerically 
Laplace inverting the frequency-dependent VACF, which 
is in turn related to the frequency-dependent friction 
through the following generalized Einstein relation given 

by 



C v (z) 



k B T 



m s (z + C(z)) 



(13) 



Thus in this scheme the frequency-dependent friction has 
been calculated self-consistently with the MSD. The de- 
tails of implementing this self-consistent scheme is given 
elsewhere. 12 ' 23 

We have evaluated the diffusion coefficient D 2 by using 
the above mentioned self-consistent scheme. The calcu- 
lated diffusion value was found to be higher than the 
simulated one. This may be partly due to the observed 
faster decay of calculated F s (k, t) than the simulated one. 
This in turn could be due to the Gaussian approximation 
for F s (k,t) which truncates the cumulant expression of 
F s (k,t) beyond the quadratic (fc 2 ) term. 29 However, the 
higher order terms which are the systematic corrections 
to the Gaussian forms can be increasingly important at 
intermediate times and wave numbers (fc). 28 

Therefore, we have performed MCT calculations using 
the simulated F s (k,t) evaluated at different values of the 
wave number, fc — nk m i ni where fc m m = 2ir/L (L stands 
for the average size of the simulation cell) and n is an 
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integer varied in the range, 1 < n < 35. F s (k,t) is then 
obtained by interpolating the simulated F s (k,t) so ob- 
tained at different wave numbers. The calculated value 
of the binary term and the density term contribution to 
the friction for different values of Mr are presented in 
table III. 

Earlier mode coupling theoretical calculations 11 for 
small solute was performed by keeping the solute and 
the solvent mass equal. In those calculations it was found 
that for size ratio Sr = 5, the solute motion is primar- 
ily determined by the binary collision between the solute 
and the solvent particles. It was also shown that due to 
the decoupling of the solute motion from the structural 
relaxation of the solvent, the contribution of the density 
mode of the solvent was much smaller than that of the 
binary term. In the present calculation we find that due 
to smaller mass of the solute, the contribution of the bi- 
nary term decreases with increase in the solvent mass as 
is clearly evident from table III. 

When compared to the simulated diffusion values, it is 
clearly seen from table II that although the MCT quali- 
tatively predicts the diffusion value for M R — 5, it breaks 
down at large values of the mass ratios. This may be be- 
cause MCT overestimates the friction contribution from 
the density mode. This breakdown of the MCT for large 
mass ratio can be connected to its breakdown observed 
near the glass transition temperature. For large solvent 
mass, the system is almost frozen and the dynamic struc- 
ture factor of the solvent decays in a much longer time 
scale when compared to the solute. So from the point of 
view of the solute, it probes an almost quenched system 
which can be expected to show the behavior very simi- 
lar to a system, near the glass transition. Just as near 
the glass transition temperature, the hopping mode also 
plays the dominant role in the diffusion process. 

V. EFFECT OF DYNAMIC HETEROGENEITY 
ON F S (K,T) OF THE SOLUTE 



good description of the dynamics of the system. The 
higher order terms which are the systematic corrections 
to the Gaussian approximation in the cumulant expan- 
sion are found to be small. 

However, in a supercooled liquid, this is not the case. 
The dynamical heterogeneities observed in a deeply su- 
percooled liquid are often manifested as the magnitude 
of the deviation of o 2 (i), the so-called non-Gaussian pa- 
rameter, from zero. It has been observed that 02(2) devi- 
ates more and more strongly and decays more and more 
slowly with increase in the degree of supercooling. 30 In 
the present study, a similar behavior has also been ob- 
served in 0:2(2) (calculated for the solute). The height of 
the maximum in 02 (t) increases as the mass of the solvent 
is increased (see figure 5), which is clear evidence that the 
solute probes increasingly heterogeneous dynamics. 

It is generally believed that the dominant corrections 
to the Gaussian result are provided by the term con- 
taining 02(2). Thus, it would be interesting to see 
whether this term alone is sufficient to explain the ob- 
served stretching in F s (k, t) of the solute at longer time. 
In order to quantify this, we have plotted in figure 10 
the simulated F s (k,t) along with the F s (k,t) obtained 
after incorporating the lowest order correction (k 4 term) 
to the Gaussian approximations for mass ratio, Mr = 
50 at the reduced wavenumber k* ~ litjoyi- For com- 
parison, the Gaussian approximation to F s (k,t) is also 
shown. It is clearly seen that the first non-Gaussian cor- 
rection to F s (k,t) is not sufficient to describe the long 
time stretching predicted by the simulation. This clearly 
indicate that at the length scales probed by the solute, 
the higher order corrections cannot be neglected. It is the 
nearly quenched inhomogeneity probed by the solute par- 
ticles over small length scales which play an important 
role in the dynamics of the system. The stretching of 
F s (k,t) observed in simulation could be intimately con- 
nected with this nearly quenched inhomogeneity probed 
by the solute particles. 



It is well-known that the self-intermediate scattering 
function, F s (k,t) can be formally expressed by the fol- 
lowing cumulant expansion in powers of fc 2 , 29 



F s (k,t) = cxp(-ifc 2 < Ar 2 (t) >) 



l + ^*2(t) 



x (ifc 2 < Ar 2 (t) >) 2 + 0(k 6 ) 



where a 2 (t) is defined as 



3 < Ar 4 (t) > 
aa( ' ) = 5<Ar 2 (t) >-- 1 



(14) 



(15) 



In the stable fluid range, it has been generally found that 
the Gaussian approximation to F s (k, t), the leading term 
in the above cumulant expansion, provide a reasonably 



VI. CONCLUSIONS 

Let us first summarize the main results of this study. 
We have investigated by using the molecular dynam- 
ics simulation the diffusion of small light particles in a 
solvent composed of larger massive particles for a fixed 
solvent-to-solute size ratio (Sr = 5) but with a large vari- 
ation in mass ratio (where the mass of the solute is kept 
constant). In addition, a mode-coupling theory (MCT) 
analysis of diffusion is also presented. It is found that the 
solute dynamics remain surprisingly coupled to the sol- 
vent dynamics even in the limit of highly massive solvent. 
Most interestingly, with increase in mass ratio, the self- 
intermediate scattering function of the solute develops a 
stretching at long time which, for intermediate values of 
mass ratio, could be fitted to a single stretched expo- 
nential function with the stretching exponent, (3 ~ 0.6. 
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In the limit of very large mass ratio, the existence of 
two stretched exponential separated by a power law type 
plateau is observed. This behavior is found to arise from 
increasingly heterogeneous environment probed by the 
solute particle as one increases the mass of the solvent 
particles. The MCT calculation of self-diffusion is found 
to agree qualitatively with the simulation results for small 
mass ratio. However, it fails to describe the simulated 
prediction at large mass ratios. The velocity correlation 
function of the solute shown interesting oscillatory struc- 
ture. 

Several of the results observed here are reminiscent of 
the relaxation of the self-intermediate scattering func- 
tion, F s (k,t) observed in the deeply supercooled liq- 
uid near its glass transition temperature. In that case 
also, one often observes combination of power-law and 
stretched exponential in the decay of the intermediate 
scattering function. We find that even the breakdown 
of MCT at large mass ratio could be connected to its 
breakdown near the glass transition temperature because 
it is the neglect of the spatial hopping mode of particles 
which is responsible for the breakdown of MCT. It is to 
be noted that those hoppings which are mostly ballistic 
in nature (after a binary collision) have already been in- 
corporated in MCT. However, MCT does not include the 
hoppings which involve collective displacement involving 
several molecules. 15 

It should be pointed out that in the MCT calculation, 
we have neglected the contribution of the current term. 
While the current contribution may improve the agree- 
ment between the simulation and MCT result for small 
mass ratio {Mr — 5), its contribution at larger mass 
ratio is not expected to change the results significantly, 
because the discrepancy is very large. 

The origin of the power law remains to be investigated 
in more detail. Our preliminary analysis shows that this 
may be due to the separation of the time scale between 
the first weakly stretched exponential (due to the disper- 
sion in the binary-type interaction term) and the second, 
later more strongly stretched exponential (which is due 
to the coupling of the solute's motion to the density mode 
of the slow solvent) . This separation arises because these 
two motions are very different in nature. However, a 
quantitative theory of this stretching and power-law is 
not available at present. 

While the origin of the stretching of F s (k,t) can be 
at least qualitatively understood in terms of the inhomo- 
geneity experienced by the solute, the origin of hopping 
is less clear. In the supercooled liquid, hopping is found 
to be correlated with anisotropic local stress 15 which is 
unlikely in the present system which is at lower density 
and pressure. 

Finally we note that the system investigated here is 
a good candidate to understand qualitative features of 
relaxation in a large variety of systems, such as concen- 
trated solution of polysaccharide in water and also mo- 
tion of water in clay. 
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TABLE I. 

The time constants (ri and T2) and the exponents {[3\ 
and P2) obtained from the stretched exponential fit to 
the F s (k, t) at the reduced wavenumber k* ~ 2n/ai2 for 
different solvent-to-solute mass ratio (Mr). 



TABLE III. 

The contribution of the binary (£,f ) and the density mode 
(i?2i) of the friction for different solvent-to-solute mass 
ratio (Mr). 



Mr = HU- 
" m 2 


C 2 B 


-DPP 


5 


23.65 


13.95 


25 


25.4 


33.8 


50 


25.45 


50.2 


250 


25.58 


99.8 





n 


0i 


T2 


02 


5 


0.082 


0.70 






25 


0.08 


0.96 


0.49 


0.67 


50 


0.083 


0.94 


0.59 


0.64 


250 


0.085 


0.91 


1.01 


0.635 



TABLE II. 



The self-diffusion coefficient values of the solute parti- 
cle predicted by the simulation and obtained from the 
MCT calculations for different solvent-to-solute mass ra- 
tio {M R ). 





r\sim 
U 2 


n MC*T 


5 


0.135 


0.1065 


25 


0.108 


0.0675 


50 


0.101 


0.053 


250 


0.0805 


0.032 
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Figure Captions 

Figure 1. The time dependence of the displacements 
for a solute particle at different solvent-to-solute mass 
ratio, Mr: (a) M R = 5, (b) M R = 25, (c) M R = 50, 
and (d) M R — 250. Note that the time is scaled by 
y/ma^/kBT. The time unit is equal to 2.2 ps if argon 
units are assumed. 

Figure 2. The same plot as in figure f , but for a sol- 
vent particle at different M R : (a) M R = 5, (b) M R = 25, 
(c) M R = 50, and (d) M R = 250. Note that here also 
the time is scaled by ^/mcr^/fcsT. 

Figure 3. The self-intermediate scattering function 
F 8 (k,t) for the solute particles for different mass ratio 
M Rl at reduced wavenumber k* — k<Jn ~ 2tt. The solid 
line represents for M R — 5, the dashed line for M R = 25, 
the dotted line for M R — 50, and the dot-dashed line for 
M R = 250. Note the stretching in F s (k, t) at longer time 
with increase in M R . For details, see the text. 

Figure 4. The self-intermediate scattering function 
F s (k,t) for the solvent particles for different mass ratio 
Mr, at reduced wavenumber k* ~ 2ir. The solid line 
represents for M R = 5, the dashed line for M R = 25, the 
dotted line for M R = 50, and the dot-dashed line for M R 
= 250. 

Figure 5. The behavior of non-Gaussian parameter 
a 2 (t) calculated for the solute particles at different mass 
ratio M R . The solid line, the dashed line, the dotted line, 
and the dot-dashed line are for M R = 5, 25, 50 and 250, 
respectively. 

Figure 6. The self-intermediate scattering function 
F s (k,t) for the solute particles for different mass ra- 
tio M R as in Fig. 3, but at the reduced wavenumber 
k* ~ 27r/c7i2. This is primarily the wavenumber probed 



by the solute. Note the emergence of a plateau at larger 
mass ratio. For details, see the text. 

Figure 7. The self-intermediate scattering function 
F s (k,t) for each of the ten individual solute particle 
obtained from a single MD run for mass ratio M R — 
250. They are calculated at the reduced wavenumber 

k* ~ 271-/(712. 

Figure 8. The plot of ln e l/D 2 vs ln e mi/m 2 , where 
D 2 is the self-diffusion of the solute, mi and m 2 are the 
masses of the solvent and solute particles, respectively. 
The slope of the straight line is about 0.13. This sug- 
gests a weak power-law mass dependence of the solute 
diffusion on the mass of the bigger solvent particles. 

Figure 9. The velocity autocorrelation function C v (t) 
for the solute particles at different values of mass ratio 
M R . The solid line represents for M R = 25, the dotted 
line for M R = 50, and the dashed line for M R = 250. The 
plot shows an increase in the negative dip with increase 
in mass ratio. For details, see the text. 

Figure 10. Comparison of the simulated self- 
intermediate scattering function F s (k,t) of the solute 
particles with the F s (k,t) obtained after incorporating 
the lowest order correction (k A term ) to the Gaussian 
approximation in the cumulant expansion (Eq. 14). The 
Gaussian approximation to F s (k,t) is also shown. The 
mean-squared displacement (< Ar(t) 2 >) and the non- 
Gaussian parameter (a 2 (t)) required as an input are ob- 
tained from the simulation. The plot is at the reduced 
wavenumber k* ~ 2-nju\ 2 and for the mass ratio M R = 
50. F s (k,t) obtained from the simulation is represented 
by the solid line, the dashed line represents the F s (k,t) 
obtained after the lowest order correction to the Gaussian 
approximation and the dotted line represents the Gaus- 
sian approximation. For the detailed discussion, see the 
text. 
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